Notebook for developing ideas to go into TellRemoval code.
Need to apply scaling of T^x to transmision of water at full resolving power and then apply a kernal to apply in at resolution of CRIRES.
Fit to the observed data (Probably with the other lines removed) to fnd the best x to apply for the correction. (Gives flatest result or zero linewidth.)
-
chip1 = "CRIRE.2012-04-07T00-08-29.976_1.nod.ms.norm.sum.wavecal.fits"
chip2 = "CRIRE.2012-04-07T00-08-29.976_2.nod.ms.norm.sum.wavecal.fits"
chip3 = "CRIRE.2012-04-07T00-08-29.976_3.nod.ms.norm.sum.wavecal.fits"
chip4 = "CRIRE.2012-04-07T00-08-29.976_4.nod.ms.norm.sum.wavecal.fits"
Obs1 = fits.getdata(chip1)
hdr1 = fits.getheader(chip1)
Obs2 = fits.getdata(chip2)
hdr2 = fits.getheader(chip2)
Obs3 = fits.getdata(chip3)
hdr3 = fits.getheader(chip3)
Obs4 = fits.getdata(chip4)
hdr4 = fits.getheader(chip4)
print("Column names = {}".format(Obs1.columns.names))
wl1 = Obs1["Wavelength"]
I1_uncorr = Obs1["Extracted_DRACS"]
wl2 = Obs2["Wavelength"]
I2_uncorr = Obs2["Extracted_DRACS"]
wl3 = Obs3["Wavelength"]
I3_uncorr = Obs3["Extracted_DRACS"]
wl4 = Obs4["Wavelength"]
I4_uncorr = Obs4["Extracted_DRACS"]
start_airmass = hdr1["HIERARCH ESO TEL AIRM START"]
end_airmass = hdr1["HIERARCH ESO TEL AIRM END"]
obs_airmass = (start_airmass + end_airmass) / 2
print("Data from Detectors is now loaded")
## Rough berv correction until input calibrated file is calibrated with non berv tapas
wl1 = wl1-.5 #including rough berv correction
wl2 = wl2-.54 #including rough berv correction
wl3 = wl3-.55 #including rough berv correction
wl4 = wl4-.7
import Obtain_Telluric as obt
tapas_all = "tapas_2012-04-07T00-24-03_ReqId_10_R-50000_sratio-10_barydone-NO.ipac"
tapas_h20 = "tapas_2012-04-07T00-24-03_ReqId_12_No_Ifunction_barydone-NO.ipac"
tapas_not_h20 = "tapas_2012-04-07T00-24-03_ReqId_18_R-50000_sratio-10_barydone-NO.ipac"
tapas_all_data, tapas_all_hdr = obt.load_telluric("", tapas_all)
tapas_all_airmass = float(tapas_all_hdr["airmass"])
print("Telluric Airmass ", tapas_all_airmass)
tapas_all_respower = int(float((tapas_all_hdr["respower"])))
print("Telluric Resolution Power =", tapas_all_respower)
tapas_h20_data, tapas_h20_hdr = obt.load_telluric("", tapas_h20)
tapas_h20_airmass = float(tapas_h20_hdr["airmass"])
print("Telluric Airmass ", tapas_h20_airmass)
try:
tapas_h20_respower = int(float((tapas_h20_hdr["respower"])))
except:
tapas_h20_respower = "Nan"
print("Telluric Resolution Power =", tapas_h20_respower)
tapas_not_h20_data, tapas_not_h20_hdr = obt.load_telluric("", tapas_not_h20)
tapas_not_h20_airmass = float(tapas_not_h20_hdr["airmass"])
print("Telluric Airmass ", tapas_not_h20_airmass)
tapas_not_h20_respower = int(float((tapas_not_h20_hdr["respower"])))
print("Telluric Resolution Power =", tapas_not_h20_respower)
#print(tapas_all_hdr)
Including the 3 tapas models to show they align well and are consistent.
plt.plot(wl1, I1_uncorr, 'b') #including rough berv correction
plt.plot(wl2, I2_uncorr, 'b') #including rough berv correction
plt.plot(wl3, I3_uncorr, 'b') #including rough berv correction
plt.plot(wl4, I4_uncorr, 'b') #including rough berv correction
plt.plot(tapas_h20_data[0], tapas_h20_data[1], "--k", label="H20")
plt.plot(tapas_all_data[0], tapas_all_data[1], "-r", label="all")
plt.plot(tapas_not_h20_data[0], tapas_not_h20_data[1], "-.g", label="Not H20")
#plt.legend()
# Make it interactive with Bokeh
bokeh.plotting.show(bokeh.mpl.to_bokeh())
(Use telluric removal modules) And plot the result.
from TellRemoval import divide_spectra, airmass_scaling, telluric_correct, match_wl
def correction(wl_obs, spec_obs, wl_tell, spec_tell, obs_airmass, tell_airmass, kind="linear", method="scipy"):
interped_tell = match_wl(wl_tell, spec_tell, wl_obs)
scaled_tell = airmass_scaling(interped_tell, tell_airmass, obs_airmass)
corr_spec = divide_spectra(spec_obs, scaled_tell) # Divide by scaled telluric spectra
return corr_spec
#def telluric_correct(wl_obs, spec_obs, wl_tell, spec_tell, obs_airmass, tell_airmass, kind="linear", method="scipy"):
# return Corrections, Correction_tells, Correction_Bs, Correction_labels
# Corrections, Correction_tells, Correction_Bs, Correction_labels = telluric_correct(wl2, I2_uncorr, tapas_not_h20[0], tapas_not_h20[1], obs_airmass, tapas_not_h20_airmass)
# Getting zero division error from this function so will try it again from te functions here
tell_airmass = tapas_not_h20_airmass
I1_not_h20_corr = correction(wl1, I1_uncorr, tapas_not_h20_data[0], tapas_not_h20_data[1], obs_airmass, tapas_not_h20_airmass, kind="linear", method="scipy")
I2_not_h20_corr = correction(wl2, I2_uncorr, tapas_not_h20_data[0], tapas_not_h20_data[1], obs_airmass, tapas_not_h20_airmass, kind="linear", method="scipy")
I3_not_h20_corr = correction(wl3, I3_uncorr, tapas_not_h20_data[0], tapas_not_h20_data[1], obs_airmass, tapas_not_h20_airmass, kind="linear", method="scipy")
I4_not_h20_corr = correction(wl4, I4_uncorr, tapas_not_h20_data[0], tapas_not_h20_data[1], obs_airmass, tapas_not_h20_airmass, kind="linear", method="scipy")
plt.plot(wl1, I1_uncorr, "b")
plt.plot(wl2, I2_uncorr, "b")
plt.plot(wl3, I3_uncorr, "b")
plt.plot(wl4, I4_uncorr, "b")
plt.plot(wl1, I1_not_h20_corr, "r")
plt.plot(wl2, I2_not_h20_corr, "r")
plt.plot(wl3, I3_not_h20_corr, "r")
plt.plot(wl4, I4_not_h20_corr, "r")
plt.plot(tapas_not_h20_data[0], tapas_not_h20_data[1], "k")
# Make it interactive with Bokeh
bokeh.plotting.show(bokeh.mpl.to_bokeh())
To use inside fit
## USEFUL functions from pedros code:
def wav_selector(wav, flux, wav_min, wav_max):
"""
function that returns wavelength and flux withn a giving range
"""
wav_sel = np.array([value for value in wav if(wav_min < value < wav_max)], dtype="float64")
flux_sel = np.array([value[1] for value in zip(wav,flux) if(wav_min < value[0] < wav_max)], dtype="float64")
return [wav_sel, flux_sel]
def chip_selector(wav, flux, chip):
chip = str(chip)
if(chip in ["ALL", "all", "","0"]):
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT1"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END4"]) # Wavelength end on detector [nm]
#return [wav, flux]
elif(chip == "1"):
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT1"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END1"]) # Wavelength end on detector [nm]
elif(chip == "2"):
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT2"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END2"]) # Wavelength end on detector [nm]
elif(chip == "3"):
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT3"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END3"]) # Wavelength end on detector [nm]
elif(chip == "4"):
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT4"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END4"]) # Wavelength end on detector [nm]
else:
print("Unrecognized chip tag.")
exit()
#select values form the chip
wav_chip, flux_chip = wav_selector(wav, flux, chipmin, chipmax)
return [wav_chip, flux_chip]
def unitary_Gauss(x, center, FWHM):
"""
Gaussian_function of area=1
p[0] = A;
p[1] = mean;
p[2] = FWHM;
"""
sigma = np.abs(FWHM) /( 2 * np.sqrt(2 * np.log(2)) );
Amp = 1.0 / (sigma*np.sqrt(2*np.pi))
tau = -((x - center)**2) / (2*(sigma**2))
result = Amp * np.exp( tau );
return result
This function broadens a spectrum assuming a Gaussian kernel. The width of the kernel is determined by the resolution. In particular, the function will determine the mean wavelength and set the Full Width at Half Maximum (FWHM) of the Gaussian to (mean wavelength)/resolution.
Parameters :
wvl : array
The wavelength
flux : array
The spectrum
resolution : int
The spectral resolution.
edgeHandling : string, {None, “firstlast”}, optional
Determines the way edges will be handled. If None, nothing will be done about it. If set to “firstlast”, the spectrum will be extended by using the first and last value at the start or end. Note that this is not necessarily appropriate. The default is None.
fullout : boolean, optional
If True, also the FWHM of the Gaussian will be returned.
maxsig : float, optional
The extent of the broadening kernel in terms of standrad deviations. By default, the Gaussian broadening kernel will be extended over the entire given spectrum, which can cause slow evaluation in the case of large spectra. A reasonable choice could, e.g., be five.
Returns :
Broadened spectrum : array
The input spectrum convolved with a Gaussian kernel.
FWHM : float, optional
The Full Width at Half Maximum (FWHM) of the used Gaussian kernel.# Convolution from Pyastronomy
from PyAstronomy import pyasl
#tapas_h20_data[0], tapas_h20_data[1], "0", 50000, FWHM_lim=5.0, plot=True
#ans = pyasl.instrBroadGaussFast(wvl, flux, resolution, edgeHandling=None, fullout=False, maxsig=None)
# Just to my CRIRES range
R = 50000
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT1"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END4"]) # Wavelength end on detector [nm]
wav_chip, flux_chip = wav_selector(tapas_h20_data[0], tapas_h20_data[1], chipmin, chipmax)
FWHM_min = wav_chip[0]/R #FWHM at the extremes of vector
FWHM_max = wav_chip[-1]/R
print("Min FWHM", FWHM_min)
print("Max FWHM", FWHM_max)
# pyasl needs equidistant wavelenghts
from Find_gcdt import gcdt
greatest_common_divisor = gcdt(wav_chip, 4)
print("Found divisor = ", greatest_common_divisor)
steps = (wav_chip[-1] - wav_chip[0])/greatest_common_divisor
print("NUmber of steps =", steps)
new_wav = np.linspace(wav_chip[0], wav_chip[-1], num=steps, endpoint=True)
# Inperolate_to_new_wave
new_flux = match_wl(wav_chip, flux_chip, new_wav)
print("length of new flux", len(new_flux))
diffs = np.diff(wav_chip)
print("First Diff={}, last diff={}".format(diffs[0], diffs[-1]))
data_step = 2 * (wav_chip[-1] - wav_chip[0])/np.min(diffs)
new_diff_wav = np.linspace(wav_chip[0], wav_chip[-1], num=data_step, endpoint=True)
new_diff_flux = match_wl(wav_chip, flux_chip, new_diff_wav)
print("length of new flux", len(new_diff_flux))
#new_flux = np.interp(new_wav, wav_chip, flux_chip)
import time
import datetime
start = time.time()
print("start time", datetime.datetime.now().time())
Conv_flux_pysal = pyasl.instrBroadGaussFast(new_wav, new_flux, 50000, edgeHandling=None, fullout=False, maxsig=None)
done = time.time()
print("end time", datetime.datetime.now().time())
elapsed = done - start
print("Convolution time = ", elapsed)
plt.plot(tapas_h20_data[0], tapas_h20_data[1],"b")
plt.plot(new_wav,Conv_flux_pysal, "r")
plt.xlabel("Wavelength (nm)")
plt.ylabel("Flux")
plt.title("Test of Pyasl convolution \n (Used interpolation to equidistant points)")
#plt.show()
# Make it interactive with Bokeh
bokeh.plotting.show(bokeh.mpl.to_bokeh())